Inter- and Intra-Chain Attractions in Solutions of Flexible Polyelectrolytes at 

Nonzero Concentration 
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Constant temperature molecular dynamics simulations 
were used to study solutions of flexible polyelectrolyte chains 
at nonzero concentrations with explicit counterions and un- 
screened coulombic interactions. Counterion condensation, 
measured via the self-diffusion coefficient of the counterions, 
is found to increase with polymer concentration, but con- 
trary to the prediction of Manning theory, the renormalized 
charge fraction on the chains decreases with increasing Bjer- 
rum length without showing any saturation. Scaling analysis 
of the radius of gyration shows that the chains are extended at 
low polymer concentrations and small Bjerrum lengths, while 
at sufficiently large Bjerrum lengths, the chains shrink to pro- 
duce compact structures with exponents smaller than a gaus- 
sian chain, suggesting the presence of attractive intrachain 
interactions. A careful study of the radial distribution func- 
tion of the center-of-mass of the polyelectrolyte chains shows 
clear evidence that effective interchain attractive interactions 
also exist in solutions of flexible polyelectrolytes, similar to 
what has been found for rodlike polyelectrolytes. Our results 
suggest that the broad maximum observed in scattering ex- 
periments is due to clustering of chains. 



I. INTRODUCTION 

Polyelectrolytes are long-chain molecules with ioniz- 
able side-groups |@,D- In a polar solvent, these ioniz- 
able groups can dissociate resulting in localized charges 
along the backbone of the polymer forming a "macroion" , 
plus mobile counterions that are no longer permanently 
bonded to the chain. Due to the long-range nature of the 
electrostatic interactions, the counterions, though mo- 
bile, are not entirely uncorrelated with the macroions. 
For polyelectrolytes with a sufficiently high charge den- 
sity along the backbone, some of the counterions remain 
bound to the chain. This phenomenon is known as "coun- 
terion condensation" . 

The qualitative theory of counterion condensation has 
been worked out by Manning || for long rodlike poly- 
electrolytes at infinite dilution. Within the Manning 
theory, counterion condensation occurs when the dis- 
tance between charges b along the chain is sufficiently 
small (i.e. the charge density is sufficiently high) com- 
pared to the length scale set by the electrostatic inter- 
actions Xb — e 2 /eksT (e being the charge, e the dielec- 
tric constant and ksT the thermal energy). Counterions 
continue to condense onto the macroions, resulting in a 
renormalized charge density along each chain, until the 



effective average charge separation b c s equals Xb- 

Manning's theory provides a lucid picture for under- 
standing counterion condensation, but its simplicity in- 
vites more complex questions. For example, it is not clear 
how the qualitative features of the theory would be mod- 
ified by having flexible chains. Equally unclear is how 
the essential results of the theory can be generalized to 
nonzero concentrations instead of just at infinite dilution. 
Recent experiments , theories |7|-[l0| and simulations 
pr| , ^2[ have attempted to address some of these issues, 
but no consistent picture has yet emerge. 

Another issue related to counterion condensation con- 
cerns the possibility of attractive interactions between 
two or more like-charged chains. There is now convinc- 
ing evidence, both numerical jl3| and analytical juj, that 
attractive interactions can exist between rodlike polyelec- 
trolytes with the same charge. This attraction is me- 
diated by the fluctuating charge density along the rods 
created by the condensed counterions when they move on 
and off the rods, leading to an effective van der Waals- 
like interaction. But the same question for flexible poly- 
electrolyte chains is still open. Since flexible chains can 
undergo conformational changes as the degree of counte- 
rion condensation changes and deviation from linearity in 
turn affects the degree of counterion condensation, these 
two factors become highly convoluted and their net ef- 
fects on the attractive interaction is unclear. 

This paper will try to address some of the issues raised 
here for solutions of flexible polyelectrolytes at nonzero 
concentrations. The complexity of the problem demands 
a careful treatment of the relevant microscopic interac- 
tions involved. Computer simulations provide a simple 
way to study the problem without resorting to approxi- 
mate methods such as the Debye-Hiickel or other mean 
field theories. In the simulations described below, all 
counterions in the solution as well as charges on the 
macroions are handled explicitly, with the full unscreened 
coulombic interactions among them. The method we 
used is similar to those in other computer simulations 
reported recently Hj. While the focus of the previ- 
ous studies has been on chain structure and single-chain 
properties, the present study concentrates on inter- and 
intra-chain interactions. 



II. MODEL AND SIMULATION METHOD 

In our simulations, each polyelectrolyte chain is repre- 
sented by a bead-spring model with N monomers. Every 
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monomer has unit charge Z = 1. We employ cubic peri- 
odic boundary condition, with M polymer chains inside 
the simulation box of dimension L 3 . In this study, we 
have considered only solutions with monovalent counte- 
rions of charge Z = — 1 and no added salt; therefore, the 
number of counterions N c necessarily equals M x N, due 
to charge neutrality. No explicit solvent molecules are in- 
cluded in the simulations. The solvent is represented by 
a dielectric continuum. The particles interact with each 
other through a coulombic potential, renormalized by the 
dielectric constant e. In addition to the coulombic inter- 
action, every pair of particles has an excluded volume 
interaction, represented by a purely repulsive truncated 
and shifted Lennard-Jones (RLJ) potential 

c/ RL.t m _ f4 £ [(f) 12 -(f) 6 + i], r<2V» ff; 

U [r) -\Q, r>2 1 /6 (J . W 

The monomer size is approximately a. We measure all 
lengths in units of a and all energies in units of e. The 
connectivity of the chains is maintained by the finitely 
extensible nonlinear elastic (FENE) potential p5[ : 

[/ FE »4^Mi-^), (2) 

which defines the spring potential between consecutive 
monomers along the chain. We set k = 7e/a 2 and 
Rq = 2a to ensure essentially no bond crossing. We 
assume the polyelectrolyte backbone is in a good sol- 
vent, therefore no explicit solvent-polymer interactions 
has been included. This model as well as the param- 
eters are essentially identical to those used in previous 
simulations by Stevens and Kremer |lT| . 

Special handling of the coulombic interaction is neces- 
sary in a solution with finite concentration. The min- 
imum image convention is usually sufficient for short 
range interactions, but the long range nature of the 
coulomb potential causes a pair of charges to interact far 
beyond their first periodic image. The Ewald summation 
method properly accounts for such contributions from all 
the images, but its implementation generally requires ex- 
tensive computation. An efficient alternative is the tab- 
ulated Ewald method with interpolation |l6j which we 
have used in the present study. This allows us to treat 
arbitrarily high concentrations with high efficiency and 
without approximation. 

The scale of the charge-charge interaction is set by the 
Bjerrum length Xb = e 2 /eksT. The simulations have all 
been carried out at a constant temperature ksT = 1.2e. 
The Bjerrum length can be varied by changing the dielec- 
tric constant e. For the results reported below, Xb varies 
between 0.83 and 6.7 a. Using the fact that the Bjerrum 
length in water at 25 °C is about 7.1 A, the typical range 
of As in an aqueous solution would be approximately 1 
to 2 a if the monomer size a is taken to be roughly 4 to 
8 A. 

We employ standard molecular dynamics (MD) meth- 
ods |l7j] for the simulation. To carry out constant tem- 
perature dynamics, either Brownian dynamics |Q or 



stochastic collision |19[ can be used. We have found 
stochastic collision to be more efficient, because it gen- 
erally generates a faster relaxation time for the polymer 
configurations |20|| . In our simulations, we assign new 
velocities to all the particles from a Maxwell-Boltzmann 
distribution every 100 to 1000 MD time steps. 

The dynamics of the system is performed using the 
Verlet leap-frog/central difference algorithm at a time 
step of 0.01 5t where r = a ^m/A8e. To ensure full equi- 
libration, the simulations were run long enough such that 
the chains move at least six times their contour lengths. 
This required anywhere from 250,000 to 4,000,000 time 
steps for most systems studied. 

We have studied multi-chain systems at three different 
monomer number densities: per 3 = 10~ 3 , 10 -4 and 10~ 5 , 
with various chain lengths N = 16, 24, 32, and 48. The 
system sizes are chosen so that the contour length of each 
chain does not exceed half of the box length. This trans- 
lates to 12 chains for AT=16, 10 for N=24, 6 or 9 chains 
for A=32, and 6 for N — 48. These systems are all well 
below the overlap density. We have also performed ex- 
tensive simulations with 27 chains at the above densities 
with N = 10. 



III. RESULTS AND DISCUSSIONS 
A. Counterion Condensation 

Counterion condensation is difficult to quantify di- 
rectly. Counting the number of counterions within a cer- 
tain distance from the polyelectrolyte is a possible though 
not very precise measure since the macroion can have in- 
fluence on a counterion far away from the chain due to 
the long-range nature of the coulomb potential. As an al- 
ternative, the osmotic pressure has often been used, both 
experimentally and theoretically, as an indirect measure 
of the degree of counterion condensation. To truly deter- 
mine the degree of counterion condensation though, it is 
necessary to study the correlation between the motions 
of the counterions and that of the macroion. This is best 
done by examining the velocity correlation between the 
counterions and the chain. 

In lieu of measuring the velocity cross-correlation, 
there is a simpler alternative for quantifying counterion 
condensation. We can measure the diffusion coefficient 
of the counterions instead. If the chains are massive, 
they diffuse slowly. Therefore, if the motions of the 
condensed counterions are correlated to the chains, they 
would exhibit a substantially reduced diffusion coefficient 
too. Measuring the diffusion coefficient of the counteri- 
ons Dc is in principle easy, at least theoretically, and 
this provides an accurate estimate of the degree of coun- 
terion condensation. In the limit where the uncondensed 
counterions are essentially ideal, i.e. assuming that they 
do not interact with each other strongly, Dc should be 
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simply proportional to the number of uncondensed coun- 
terions. The counterion diffusion coefficient Dq can also 
be related to the mobility pic of the counterions via the 
Einstein relationship. 

The counterion diffusion coefficients Dc are shown in 
Fig. |] for various Bjerrum lengths at the three different 
concentrations studied, closed triangles for pa 3 = 10~ 5 , 
open squares for 10 -4 and closed circles for 10~ 3 . In 
order to estimate the intrinsic mobility of the counteri- 
ons in the absence of counterion condensation, we have 
also run the simulations at these three concentrations 
for various Xb without the macroions. The dashed line 
in Fig. [l] indicates this intrinsic diffusion coefficient D c , 
which was found to be independent of A^ for this concen- 
tration regime and roughly independent of p also. The 
open circle at A^ =0 indicates the diffusion coefficient 
obtained for a system with the coulombic interactions 
turned off, providing another independent estimate for 
the intrinsic diffusion coefficient D c in the absence of 
counterion condensation. 
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FIG. 1. Self-diffusion coefficient Dc of counterions at 
three different concentrations pa 3 = 10 -3 (closed circies), 
10 -4 (open squares) and 10 -5 (closed triangies) as a func- 
tion of Bjerrum length As. For an estimate of the intrin- 
sic self-diffusion coefficient in the absence of any counterion 
condensation, the open circie indicates Dc obtained from a 
simulation with the coulombic interactions turned off and the 
dashed line indicates Dc from simulations with counterions 
only but no polyelectrolyte chains. 

As expected, Dc decreases with increasing Xb, sig- 
naling an increasing counterion condensation. For the 
lowest concentration studied pa 3 = 10 -5 (triangles in 
Fig. the dependence is apparently linear, and this lin- 
earity appears to extend throughout the entire range of 
As studied. For still larger values of A^ which we have 
not studied, this linearity is expected to eventually give 



way to a faster than linear decrease, since Dq obviously 
cannot be less than zero. For the two higher concentra- 
tions, the dependence deviates from linearity sooner, but 
the basic features are the same. 

There is no sharp break in the Dc dependence on As. 
This suggests a behavior qualitatively different from that 
predicted by Manning theory. Manning theory suggests 
that there is some critical value of Xb = b below which 
no condensation occurs. For the model used here, b, the 
mean distance between charges on the chain, is 1. Ac- 
cording to Manning theory, the onset of counterion con- 
densation should then occur at A^ = 1. However, we 
observe in Fig. [l] that for all concentrations, the depen- 
dence of Dc on As does not exhibit any break at As near 
1 or at any other value. 
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FIG. 2. Self-diffusion coefficient Dc of counterions in solu- 
tions containing partially-charged polyelectrolyte chains with 
charge fraction /. Three different concentrations are shown: 
pa 3 = 10 -3 (closed circies), 10 -4 (open squares) and 10 -5 
(closed triangles) and As = 1.7. 

This departure from Manning theory is likely a result 
of the flexibility in the chain conformation, which in turn 
is strongly dependent on the degree of counterion con- 
densation. As we will show in Sect. Ill B , the structure 



of flexible polyelectrolytes are not always extended, and 
this affects the validity of the arguments behind Manning 
theory, which was originally formulated for long rods. 

Using the simulations, we can examine other qual- 
itative aspects of Manning theory. For example, one 
can study the condensation of counterions onto partially- 
charged polyelectrolyte chains. The results presented in 
Fig. |l| are for fully-charged chains, i.e. every monomer 
on the backbone is charged. We have also examined 
partically-charged chains. By charging only certain frac- 
tions of the monomers, we simulated solutions of chains 
with different charge fractions /. The condensation of 
counterions was again quantified by measuring the coun- 
terion diffusion coefficients D c . According to Manning 
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theory, counterfoil condensation will continue until the 
renormalized charge fraction /„ (charge of the backbone 
minus the charge of the condensed counterions divided 
by the number of monomers) on the chain is such that 
the average charge separation is equal to the Bjerrum 
length Xb- Consequently, Manning theory predicts that 
the renormalized charge fraction /„ as a function of the 
charge fraction of the bare chain / would exhibit a sat- 
uration behavior. Figure ^ shows how D c depends on 
/ for three different concentrations at Xb = 1.7. From 
these data we can determine the number of condensed 
counterions and compute the renormalized charge frac- 
tion /„ . Figure |3| shows /„ extracted from results in Fig. ^ 
plotted as a function of /. For this set of parameters, 
the renormalized charge fraction f n should saturate at 
0.6 according to Manning theory. But the data show no 
such saturation. Instead, the renormalized charge frac- 
tion f n seems to increase monotonically for all concen- 
trations studied. 



incorrect. 
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FIG. 3. The renormalized charge fraction f n on each chain 
as a function of the charge fraction on the bare chain extracted 
from counterion diffusion coefficient data from Fig. ^. Three 
different concentrations are shown: pa :i = 1CF 3 (closed cir- 
cles), 10 -4 (open squares) and 10 -5 (closed triangles) and 
A_b = 1.7. The data do not exhibit the saturation behavior 
predicted by Manning theory. 

Our observation that there is no saturation in the 
renormalized charge fraction is consistent with the an- 
alytical theory of Liu et al. |l(| for finite-length polyelec- 
trolyte rods in a poor solvent. However, a recent small- 
angle X-ray and neutron scattering experiment |(| seems 
to suggest that saturation does occur, in support of Man- 
ning theory. This experimental conclusion was based on 
the assumption that the peak observed in the scatter- 
ing amplitude arose from a correlation tube of the size 
of the Debye-Huckel screening length around each chain. 
Since the screening length is thought to depend on the 
concentration of free counterions, it was argued that the 
position of the scattering peak could be used to deter- 
mine the renormalized charge fraction on the chains. In 
Sect. Ill C\ we will show that this interpretation may be 



1.0 



D 



0.5 



0.0 




-5.0 



-4.0 



-3.0 



log po 



FIG. 4. Data from Fig. replotted as a function of con- 
centration for different Bjerrum lengths \b = 0.83 (triangles), 
1.7 (diamonds), 3.3 (squares) and 6.7 (circles). 
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FIG. 5. The fraction of condensed counterions extraced 
from data in Fig. ^| for different Bjerrum lengths \b ~ 0.83 
(triangles), 1.7 (diamonds), 3.3 (squares) and 6.7 (circles). 

The counterion diffusion coefficient in Fig. [I] are re- 
plotted in Fig. |] as a function of concentration for the 
four different Bjerrum lengths studied. At each Xb, there 
is a steady decrease in Dc with increasing concentra- 
tion, indicating that counterion condensation increases 
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with increasing concentration. Again, the dependence 
is smooth, without any abrupt changes or any signs of 
a critical transition point. Our results are qualitatively 
consistent with the experimental measurements of coun- 
terion self-diffusion coefficients using pulsed field gradi- 
ent NMR techniques in the low concentration region [£lj . 
From these data, we can directly determine the fraction 
of counterions condensed onto the chains. This is shown 
in Fig. U as a function of concentration for the four dif- 
ferent Xb studied. Clearly, the number of counterions 
condensed is not only a function of the Bjerrum length, 
but also a strong function of concentration. 
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FIG. 6. Log-log plot of the average squared radii of gyra- 
tion as a function of chain length TV for different Bjerrum 
lengths and different concentrations pa 3 — 10~ 3 (circles), 
10~ 4 (squares) and 10 -5 (triangles). The dashed lines indi- 
cate results from simulations with the coulombic interactions 
turned off. 



B. Single-Chain Conformations and Intrachain 
Attractions 

We have studied the single-chain structure at vari- 
ous concentrations and Bjerrum lengths for a number 
of different chain lengths N. The equilibrium averaged 
squared radii of gyration (R 2 ,) are shown in Fig. ^| for 
four different Bjerrum lengths. Data for the three differ- 
ent concentrations are displayed, with triangles for pa 3 
= 1CP 5 , squares for ICP 4 and circles for 10~ 3 . For all 
three concentrations, the data fall on straight lines on a 
log-log plot for all Xb studied, indicating that the chain 
lengths from the shortest one with N = 16 on are all in 
the scaling limit. For comparison, we have also shown 
the squared radii of gyration for neutral chains at infi- 
nite dilution as dashed lines in Fig. [| For all lengths 
studied, the charged chains are substantially larger than 
the neutrals. 

From the slopes of the lines in Fig. || we can extract 
the exponent v in the scaling relationship (R 2 g ) ~ N 2u , 
and their values are shown in Table Q. For the lowest 
concentration pa 3 — 10~ 5 , the chains are fully extended 
with v w 1 for all Xb ■ For the intermediate concentration 
pa 3 = ICP 4 , the chains are still fairly extended, with a 
slight variation in v going toward larger A^ . 

For the highest concentration pa 3 = 10 -3 , however, 
the variation in the exponent v is quite substantial. For 
the smallest Xb = 0.83, v w 0.85, indicating a semi- 
extended chain structure. But at the largest A^ = 6.7, 
v w 0.48, indicating a globular structure with a scaling 
exponent that is actually smaller than a gaussian chain. 
In between these two extremes, there is a consistent de- 
crease in v, suggesting that the size of the chain decreases 
steadily toward more globular structures. Figure shows 
examples of chain conformations at pa 3 = 10~ 3 for dif- 
ferent Bjerrum lengths, illustrating the progression from 
a semi-extended structure to a globular structure in ac- 
cordance with the variation in v. Though the scaling 
indicates an exponent smaller than gaussian behavior for 

TABLE I. Scaling exponents for (-Rg) at various Bjerrum 
lengths and concentrations. 



Xb p v 

0.833 W 3 084 

10" 4 0.98 
10~ 5 1.0 

1.667 10" 3 0.77 

10~ 4 0.99 
10 -5 1.0 

3.333 10" 3 0.62 

10 -4 0.94 
10~ 5 1.0 

6.667 10~ 3 0.48 

10~ 4 0.82 
10 -5 1.0 
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large values of Xb , it is important to remember that the 
absolute size of the charged chains are larger than the 
neutrals (dashed lines in Fig. O). 
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FIG. 7. Typical chain conformations for concentration per 3 
= 10 -3 at four different Bjerrum lengths, As: (a) 0.83u, (b) 
1.7cr, (c) 3.3cr and (d) 6.7<r. Light-colored spheres represent 
monomer units on the polyelectrolyte chains and dark-colored 
spheres the counterions. (Differences in the size of the spheres 
are for display purpose only. The counterion and monomers 
have identical effective size in the simulations.) Notice the 
variations in the degree of counterion condensation as a func- 
tion of Xb ■ 

Recent studies of effective attractive interactions be- 
tween like-charged polyelectrolytes have considered only 
rigid rods |l3| , |l4| ] where the issue of intrachain attraction 
did not come up. A natural extension of these studies 
is to ask whether an effective attractive interaction also 
exists between different parts on the same chain, and if 
it docs, whether its nature is identical to the interchain 
interactions. The globular structures shown in Fig. |?j 
for large Bjerrum lengths suggest the existence of at- 
tractive intrachain interactions. This is corroborated by 
the smaller-than-gaussian scaling exponent v observed. 
Compared to the extended structures at small Xb, the 
higher Xb globular structures are more contorted, with 
a larger number of bends. The same behavior has also 
been observed in a previous simulation of flexible poly- 
electrolytes jll] . The explanation ascribed to it was that 
increasing counterion condensation at larger Xb screened 
the intrachain repulsion in such a way that the structure 
of the chains went from being fully extended to that of 
an essentially neutral polymer. But the scaling evidence 
from the present work goes beyond this assertion and 



suggests the presence of an attractive interaction caus- 
ing the chains to attain a structure even more compact 
than that of a neutral polymer. This collapse of the chain 
agrees qualitatively with the recent results of Winkler et 
al. @ 

Having established the presence of an attractive in- 
trachain interaction, the next question is whether this 
attractive interaction has the same origin as the inter- 
chain attraction observed in polyelectrolyte rods. This is 
a difficult question to answer. Interchain attraction can 
arise from the charge fluctuations along the backbone of 
the rods as condensed counterions come on and off the 
rods. But this is not the only mechanism possible. Ray 
and Manning recently suggested that two like-charged 
rods can share counterions in such a way to form what 
is analogous to a covalent bond p^| . The effective dis- 
tance scale of this covalent-like interaction is predicted 
to be of the order of the Debye-Huckel screening length, 
which is much longer than the intrachain distance here. 
Therefore, the first mechanism is more likely to be cor- 
rect for the attractive intrachain interactions. However, 
more work is necessary to fully elucidate the origin of the 
intrachain attraction. 

The attractive interactions are evident only for solu- 
tions with very large Bjerrum lengths. For smaller Xb, 
the chains are apparently well extended. For these cases, 
the intrachain attractions, if present, are not strong 
enough to overcome the coulombic repulsion between dif- 
ferent parts of the same chain. 



C. Interchain Interactions 

The most straightforward way to determine whether 
interchains attraction exists is to calculate the free energy 
as a function of the separation between two chains. The- 
oretically, this could be done by an umbrella sampling 
procedure |23| or by computing the potential of mean 
force of separating two chains |24j| . Unfortunately, due 
the long-range nature of the coulomb potential, both of 
these methods suffer from heavy statistical noise. Com- 
pounded by the need to deal with a finite polymer con- 
centration (i.e. multiple chains), these direct methods 
are not always very useful in practice. 

An alternative is to monitor the radial distribution 
function g(r) of the centers-of-mass of the polyelectrolyte 
chains. The logarithm of g(r) is related to the free en- 
ergy of separating a pair of chains, so by monitoring g(r) 
we can obtain information equivalent to those given by 
direct calculations of the free energy. From sufficiently 
long MD simulations carried out at equilibrium, an ac- 
curate determination of the radial distribution function 
is not too difficult. 

The data we will show were obtained from very long 
simulations of M = 27 chains with N = 10 monomers 
at many different concentrations and Bjerrum lengths. 
These chains are shorter than those used for obtaining 
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the results shown in the last two subsections. The moti- 
vation for this was to enhance the statistics in the g(r) 
calculations by putting more chains in the simulations, 
but reducing the size of each chain to keep the computa- 
tion at a manageable level. We have also run calculations 
at a number of different combinations of M and larger N. 
The results are all qualitatively identical to those shown 
below for M = 27 and N = 10. 

Figure || shows the radial distribution function g(r) 
as a function of the center-of-mass separation between 
chains for a very small Bjerrum length A^ = 0.83cr. The 
simulations were carried out for three different concen- 
trations by varying the size of the simulation box L: pa 3 
= 1CT 5 (L/a = 300), 10~ 4 (L/a = 139) and 10~ 3 (L/a 
= 64). 
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FIG. 8. The polymer radial distribution function for the 
centers-of-mass of the chains at Bjerrum length Xb — 0.83a, 
for three different concentrations pa 3 — 1CP 5 (L/a = 300), 
10" 4 (L/a = 139) and 10" 3 (L/a = 64). 

Although we have used very long runs to extract the 
g(r) functions, there is still substantial statistical noise 
in the data. We have indicated the size of the statistical 
noise by error bars which represent one standard devi- 
ation from the mean. Despite this, some of the minor 
peaks and valleys in g(r) may not be significant. So we 
will rely on the gross features of g(r) to draw our conclu- 
sions. 

For the lowest concentration (L/a — 300), g(r) is fairly 
featureless with two broad humps. There is a correlation 
hole for distances smaller than roughly lOOer, indicat- 
ing that chains cannot come closer together by distances 
much smaller than lOOer. Beyond that, there is a fairly 
constant density. Apparently, for this A^ at this concen- 
tration, the polyelectrolyte chains behave very much like 



particles with an effective radius of approximate lOOer. 
Because there are 27 chains in the periodic box of length 
L = 300cr, lOOer happens to be the dimension of the av- 
erage volume available for each chain. So the effective 
radius of lOOer deduced from the g(r) appears to be a 
result of a net repulsive interaction pushing chains away 
from each other to fully fill the volume available. It is 
important to point out here that the radius of gyration 
of the chains are much shorter than the dimension of 
the correlation hole and the concentrations studied here 
are all well below the overlap density. Incidentally, the 
Debye-Hiickel screening length for these conditions is also 
roughly 100u. Therefore, it is also possible that the corre- 
lation hole in g(r) may be related to the screening length 
instead. That this is actually not the case will be evident 
from data to be presented below. 

Now turning to the intermediate concentration pa 3 = 
10~ 4 (L/a =139), we observe similar behavior as for the 
lower concentration. This time, the dimension of the cor- 
relation hole shrinks to about 50a, which again is about 
1/3 of L. Beyond this, g(r) is again rather featureless, 
indicating an essentially structureless solution of chains 
with repulsive interactions. 




FIG. 9. Same as Fig. g for \ B = 1.7a. 

Next, we examine the highest concentration pa 3 = 
10~ 3 (L/a = 64). g(r) shows a distinct peak in the 
short-distance region. Beyond that, g(r) has the gener- 
ally featureless appearance similar to the last two concen- 
trations. Based on the observations in the previous two 
cases, we expect the dimension of the correlation hole to 
be at 1/3 the box size, namely around 20a. But the left 
edge of g(r) now stands at approximately r — 10a. These 
facts together indicate that the peak at small distance is 
a new feature not present in the previous cases, and the 
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left edge of g(r) does not actually correspond to the di- 
mension of the correlation hole. The origin of this new 
feature will become clear below. 

In Fig. ||, we show the same data for a larger Xb = 
1.7a. For the lowest concentration, a similar g(r) is ob- 
served here compared to the one at the smaller Bjerrum 
length in Fig. ||. Notice that the dimension of the correla- 
tion hole has not changed, although the screening length 
is smaller for Xb — 1.7a than it is for 0.83er in Fig. || by a 
factor of \[2. So the size of the correlation hole is appar- 
ently unrelated to the Dcbye-Hiickel screening length. 

The size of the correlation hole appears to be more 
closely related to the polymer concentration instead of 
the screening length. Essentially the same conclusion 
can be inferred from the study of Stevens and Kremer 
jnj, which showed that in the very low concentration 
regime, the position of the peak in the structure factor 
scales with the polymer concentration to a power of ap- 
proximately 1/3. In this very dilute regime, the peak 
in the structure factor is entirely due to the correlation 
hole in g(r); therefore, the size of the correlation hole 
should scale with the polymer concentration to a power 
of —1/3. This is corroborated by our observations here 
that the correlation hole shrinks by a factor of roughly 2 
when the density decreases by a factor of 10. 

The fact that the size of the correlation hole has no 
apparent relationship to the screening length has impor- 
tant implications for the small-angle X-ray a nd ne utron 
scattering experiments || referred to in Sect. Ill A. Our 
results suggest that the arguments used in the experi- 
ments in arriving at the conclusion that the counterion 
condensation saturates at a certain renormalized charge 
fraction is wrong, i.e. the size of the correlation hole 
extracted from the peak in the scattering intensity may 
have nothing to do with the screening length at all. 

For the intermediate concentration in Fig. ^, there ap- 
pears to be an extra feature for distances smaller than the 
expected size of the correlation hole (50c as in Fig. ||) just 
like in the highest concentration in Fig. ||. Given the size 
of the statistical noise and the quality of the data, we are 
uncertain about the significance of this feature. But the 
same feature will become more prominent in the data to 
be shown below. So we skip to the highest concentration. 

For the highest concentration in Fig. ^, the feature 
observed previously for Xb = 0.83a is now much more 
pronounced. There is now a sharp peak at r w 10a. 
This peak is indicative of some kind of clustering of two or 
more chains at short distances. By integrating under this 
peak, we can estimate the average number of neighbors 
each chain has. This turns out to be approximately 0.41. 
This is strong evidence that there exists an attractive 
interaction between chains, causing them to cluster with 
each other. Indeed, if we examine the configurations, 
we observe a small number of chains pairing up. These 
pairs appear to be true equilibrium structures, because 
they would dissolve and reform continually, even in cases 
where the short-distance peak in g(r) cannot be clearly 
resolved. For example, the integrated peak height in the 



radial distribution function up to r 
10 for pa 3 



Fig. 



= 10" 



22.5a is plotted in 
4 and Xb — 0.83er, indicating that 



the clustering is a transient but equilibrium phenomena. 



0.40 



0.20 



0.00 




0.0 500.0 1000.0 1500.0 

thousand time steps 

FIG. 10. Equilibrium fluctuations of the integrated num- 
ber of chains under the first peak in gjr) up to a cutoff dis- 
tance of 22.5<r from the data of Fig. pi at pa 3 = 10 -4 as a 
function of MD time steps. 

Figures [ll] and [l2] show the same data for two larger 
Bjerrum lengths Xb — 3.3cr and 6.7a. For the two high- 
est concentrations, the sharp peak at short distance con- 
tinues to grow in intensity. The feature that was barely 
visible for the intermediate concentration in Fig. ^|is now 
fully developed. The rise in the short-distance peak in- 
dicates that more chains are clustering as either the con- 
centration or the Bjerrum length or both increase, im- 
plying that the attractive interactions grow in strength 
at the same time. The rising peak in g(r) suggests that 
more chains are participating in these clusters. Again, 
by integrating under the short-distance peak in g(r), we 

TABLE II. Average number of neighbors around each 
chain obtained from integrating under the first peak in the 
radial distribution function up to cutoff distance r c . 



p 


r c /cr 


X B 


n(r) 


10" a 


9.75 


0.867 


0.29 






1.667 


0.41 






3.333 


1.69 






6.667 


1.88 


io- 4 


25.5 


0.867 


0.32 






1.667 


0.39 






3.333 


0.74 






6.667 


0.58 
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can estimate the average number of neighbors each chain 
has, and the results are tabulated in Table [n] for the two 
highest concentrations. For instance, the peak in g(r) 
for As = 6.7 and concentration per 3 = 10~ 3 indicates 
that each chain has approximately 1.8 neighbors on the 
average within a distance of 10a around it. 




FIG. 11. Same as Fig. N for X B = 3.3a. 




FIG. 12. Same as Fig. g for X B = 6.7a. 

Notice that in spite of the sharp peak observed at short 
distance, there are no other distinct secondary peaks. If 
a large number of chains come together to form a large- 



scale structure, there would have been second and third 
neighbor peaks evident in g(r). The absence of these 
peaks suggests that most of the intensity under the sharp 
peak at short distance is in fact due to small clusters, in 
agreement with the data in Table [fi|. Indeed, a visual ex- 
amination of the configurations indicates that most of the 
chains are forming pairs or triplets. Examples of pairs are 
shown in Fig. The precise reason for the absence of 
larger clusters is not clear at this point. There are several 
possibilities: (a) Some special feature of the attractive 
interchain interaction causes chains to form small clus- 
ters but precludes the creation of larger clusters; (b) The 
number of chains used in the simulations is too small to 
accurately reflect the presence of large-scale structures; 
or (c) The lengths of our simulations are not long enough 
for the true equilibrium large-scale structures to emerge. 
Larger simulations are needed to determine the actual 
reason for the apparent absence of large clusters. 




FIG. 13. Examples of chains pairing at two different Bjer- 
rum lengths Xb = (a) 3.3a and (b) 6.7a at the highest con- 
centration pa 3 = 10 -3 . Notice the high degree of interchain 
correlations. 

Also interesting is the Xb dependence of the position 
of the small-distance peak. This is most clearly seen for 
the highest concentration. Starting from Xb = 0.83cr, the 
peak position shifts consistently to smaller values, from 
about 10a to 6a at Xb = 6.7a. A similar trend can be 
seen in data for the intermediate concentration. 

The data in the radial distribution function demon- 
strate unequivocally that attractive interchain interac- 
tions do exist in solutions of like-charged flexible polyelec- 
trolytes. The strength of this attraction also grows with 
polyclcctrolytc concentration and the Bjerrum length. 
These attractive interactions lead to clustering of chains, 
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which should be observable in neutron and light scatter- 
ing experiments. Indeed, it has been suggested that the 
broad peak detected in recent neutron and light scatter- 
ing experiments B at small wavevectors may actually be 
related to the formation of clusters or domains || . 

To verify the connection between the short-distance 
peak observed in the polymer center-of-mass radial dis- 
tribution function and the broad low-wavevector peak 
detected in the scattering intensity we can compute the 
structure factor S(q) from the simulation data and com- 
pare them to the experiments. Unfortunately, due to 
the small size of the simulation box used, the number 
of allowed wavevectors dictated by the periodic bound- 
ary condition we have used is extremely sparse at the 
low-wavevector region. Instead of Fourier transforming 
the monomer density directly to get S(q), we first com- 
puted the monomer-monomer radial distribution func- 
tion g m (r) and then S(q) from it using the formula 
S(q) = 1 + p m h(q), where h{q) is the Fourier transform 
of h(r) = g m (r) — 1 and p m is the monomer density. We 
have compared the S(q) computed this way with a direct 
Fourier transform of the density, and found the qualita- 
tive features to be identical. 



15 



I 10 

s _ 
a* 5 



the part of the structure factor S(q) that arises solely 
from the interchain correlation, and this is shown in 
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FIG. 14. 



Total monomer-monomer radial distribution 



function g m (r) for As = 1.7 and pa = 10 . Inset shows 
the interchain contribution to g m (r) (solid line) and the poly- 
mer center-of-mass radial distribution function (circles) taken 
from Fig. [| 

Using the monomer radial distribution function offers 
another advantage. It allows us to separate the inter- 
and intra-chain contributions to the total scattering in- 
tensity. Figure R4] shows the total monomer g m {f) for 
\b = 1.7 and pa — 1CP 3 . Clearly, the total (intra-chain 
plus inter-chain) g m (r) is dominated by the intrachain 
correlations at short distance. The inset shows only the 
inter-chain contribution to g m {r) at a much expanded 
scale as the solid line. Also in the inset, superimposed 
on the monomer g m {f) is the polymer center-of-mass g(r 
for the same Xb and concentration taken from Fig, 
Clearly, the interchain contribution to g m (r) closely re- 
sembles the polymer g(r) we discussed earlier. From the 
interchain contribution to gm{r), we can then compute 



Fig. 15. The S(q) shows a broad peak with a maximum 
at g max « 0.18cr -1 . By truncating the Fourier integral 
at different upper limits, we can verify directly that the 
peak in S(q) is essentially due to the first and second 
peaks in g m (r). We found similar results for all other 
concentrations and Bjerrum lengths. 
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FIG. 15. Structure factor derived from the interchain 
monomer radial distribution function in Fig. [l4] showing broad 
peak similar to those observed in scattering experiments. 

Although we have correlated the peak in the structure 
factor to the attractive short-distance peak in the poly- 
mer g(r), the accuracy of our data precludes a reliable 
determination of how <7 ma x scales with the polymer con- 
centration. The reason for this is that although we have 
seemingly circumvented the problem of the sparseness of 
the allowed wavevectors by analyzing the monomer radial 
distribution function, the periodicity in principle does not 
permit us to extract real independent information about 
S(q) for wavevectors q which do not correspond to recip- 
rocal vectors of the cubic boundary condition we employ. 
If the actual peak lies between two reciprocal vectors, 
we would not be able to accurately measure its position. 
Given the size of the simulation box we have used for 
these calculations, the reciprocal vectors are certainly too 
sparse. A much larger simulation would be required for 
a direct comparison with experimentally obtained scat- 
tering intensities. 

Finally, we can compare the general characteristics of 
the attractive interactions seen in our simulations of flex- 
ible polyelcctrolytcs with those of the rigid rods observed 
in previous theoretical work fr[3 14|. First, the attractive 



interactions between the flexibles are strongly concen- 
tration dependent whereas for the rods they exist even 
for infinite dilution. This is not really very surprising, 
given the fact that the attractive interactions are results 
of charge fluctuations due to the condensed counteri- 
ons. For infinite rods, condensation could occur even 
at zero concentration, but for flexibles it could not. The 
degree of counterion condensation for flexibles increases 
with concentration, and so should the strength of the at- 
tractive interactions. Secondly, the distance scale of the 
attractive interactions seems to be very different in the 
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flexibles compared to the rods. For rods, previous cal- 
culations show that the attractive interactions are very 
short-ranged, of the order of 10 or 20 A. For the flexibles, 
if we assume a monomer size of roughly 4 A, the typi- 
cal position of the peaks observed in g{r) from Figs. || 
to |l2j would be anywhere from 30 to 100 A, depending 
on the concentration and the Bjerrum length. On the 
other hand, these distances scales are very much in ac- 
cord with those suggested by neutron and light scattering 
experiments 

IV. CONCLUSIONS 

We have carried out detailed molecular dynamics sim- 
ulations to study inter- and intra-chain interactions in 
solutions of flexible polyelectrolytes at nonzero concen- 
trations with no added salt. The simulations were per- 
formed with explicit counterions and unscreened coulom- 
bic interactions. Specifically, we investigated questions 
related to the existence of effective attractive interactions 
as well as the nature of counterion condensation at finite 
concentration. 

From the counterion self-diffusion coefficient, we were 
able to determine how counterion condensation depends 
on both the Bjerrum length and the polymer concen- 
tration. The data revealed a condensation mechanism 
that is distinctly different from that predicted from Man- 
ning theory. Specifically, we observed no saturation in 
the renormalized charge fraction on the chains which ap- 
peared to be a smooth function of the Bjerrum length as 
well as the polymer concentration. 

From the scaling of the mean squared radii of gyration 
with the chain length, we extracted the scaling expo- 
nent under a number of different conditions with various 
Bjerrum lengths and concentrations. For large Bjerrum 
lengths and high polymer concentrations, a scaling expo- 
nent smaller than a gaussian chain was observed, indicat- 
ing that the chains had a more compact structure than 
ideal chains. In conjunction with a visual examination of 
the chain conformations, we concluded that the compact 
structures are results of attractive interactions between 
different parts of the same chain. 

To ascertain the presence of attractive interchain in- 
teractions, we have studied in great detail the radial dis- 
tribution functions of the centers-of-mass of the chains in 
solution. We found a peak at small distances in the radial 
distribution function for high Bjerrum lengths and/or 
high concentrations. From the lack of secondary peaks 
in the radial distribution function and a careful exam- 
ination of the chain configurations, we concluded that 
this peak is a result of small domains formed by chains 
clustering in the solution. This behavior is indicative of 
the presence of strong interchain attractions. From the 
structure factor, we have directly verified the connection 
of the small-distance peak in the polymer radial distribu- 
tion with the broad peak detected at small wavevectors in 



the experimental small-angle X-ray and light scattering 
intensities. 
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